#include <SeqJacobi.hpp>
#include <algorithm>
#include <cmath>

// Construtor
SeqJacobi::SeqJacobi 	(std::vector<std::vector<double> >& A,
				 		 std::vector<double>& B,
				 		 std::vector<double>& x,
				 		 double error,
				 		 int maxIterations) :
A(A), B(B), x(x), error(error), maxIterations(maxIterations), solved(false)
{
	// Nothing to be done here.
}

std::vector<double>
SeqJacobi::solve(int * iterations) {
	if (solved) return x;

	std::vector<double> vx[2];
	vx[0] = x;
	vx[1] = x;
	int curx = 0;
	int myIterations = 0;
	int n = A.size();

	// Prepara a matriz
	for (int i = 0; i < n; i++) {
		double diag = A[i][i];
		for (int j = 0; j < n; j++) 
			A[i][j] /= diag;
		B[i] /= diag;
		A[i][i] = 0.0;
	}

	// Enquanto nao chegar ao limite de iteracoes e nao convergir
	for (bool biggerThanError = true;
		 myIterations < maxIterations && biggerThanError; 
		 myIterations++) {

		// Calcula x^(k+1)
		for (int i = 0; i < n; i++) {
			vx[curx][i] = B[i];
			for (int j = 0; j < n; j++)
				vx[curx][i] -= A[i][j] * vx[!curx][j];
		}

		// Colhe valores para checar convergencia
		double maxDiff = 0.0;
		double maxX = error;
		for (int i = 0; i < n; i++) {
			maxX = std::max(maxX, std::abs(vx[curx][i]));
			maxDiff = std::max(maxDiff, std::abs(vx[curx][i] - vx[!curx][i]));
		}

		// Verifica se convergiu
		if (maxDiff / maxX <= error) biggerThanError = 0;

		// Alterna entre x^k+1 e x para efetuar nova iteracao
		curx = !curx;

	}


	if (iterations) *iterations = myIterations;

	x = vx[!curx];
	solved = true;
	return x;
}
